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APPENDIX c 

LEAD 3 COMPUTER PROGRAM 


)ISCUSSI 0N OF TH E_PROGRAM to 

a colter P^ra™ vh.ch c.a _ 

flux and flux divergence through a nonrs 
determine the radiatrve ^ ^ ^ conside rs species topical of 

thermal planar slab o g ^ carbon phenolic ablation 

alr at high temperatures an ^ ^ ^ depth llne 

a r't-c The equations solve 

products. T q ce ,rtions and line widths 

using species cross sections 
and continuum gas model using 

as the basic data. . 

Th e detailed radiation model used in the present 

, i rWelooed by Wilson (Ref. C.l). 
coupled line and continuum model develop y ^ ^ } ^ 

ter program used in the present analysis i 

^ . n o£ Wilson's subroutine TRAMS. »AB 3 consists 

ls a modified version of Mils 

and a set of subroutines which are for 
O* • **»* prOSr “ a . dr iver program and 

eart subsets of subroutine TRAMS. By providing 

set of subroutines greater flexabllity an 
breaking TRANS into a set of 

reduced computation time was achieve 

c i fnol for evaluating tn 
«- program TRAP 3 provides a useful ^ 

radiative flux and flux divergence across a s ^ 

both air and ablation species. The program can be used 

. ..h a flow- field program or independently 

form in conjunction with (r eque ncy bands 

parametric studies. Twelve continuum 

are used in the radiation calculation. Radiation properties 
following species are included. 
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I.he computational techniques and equations solved are presented 
„ detail in Appendix B. The radiation model and species considered 
,ere developed as a compromise between detail of computation (computer 
time) and accuracy. The accuracy of the line and continuum calculation 
for the four atomic species was found to be quite good by Wilson (Ref. 
U.l) by comparison with a more detailed radiation program RATRAF. 

INPUT GUIDE 

All inputs to the radiative transfer computer program (LRAB 3) are 
read from cards-, no tapes are required. The basic inputs consist of 
(1) temperature: (2) number densities or mole fractions, density and 

average molecular weight: (3) shock layer coordinates: (4) input 

option parameters-, and (5) a set of shock layer points used in the 
integration to obtain the flux. Three basic formats are used for 
input, 15, E12.0 and A4. The 15 is an integer format consisting of 
five (5) columns right adjusted, the E12.0 is a floating pornt format 
occupying twelve (12) columns with a decimal point £unched on the 
card, and A4 is an alphabetic format. Multiple cases may be run by 
placing the input data for each new case behind the data for the 

previous one. 

Tab. C.l provides the details of the card input and Tab. C.2 
gives the meanings of the input variables. 
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CARD TYPE 

1 

2 

3 

4 

5 


TABLE C.l 

CARD INPUT FOR LRAD 3 
Variables 

TIT (I) 

NETA, MF, NS, LINES, IDG, IEZ 
R, DELTA, DTIL, XM0L, RDZ 

T (I) 

YD (I) 


6 

ETA (I) 

7 

INDEX 

8A* 

DENS (K, INDEX) 

8B* 

FRAC (K, INDEX) 

9° 

ETZ (I) 

lo- 

RH0 (I) 

ll* - 

AMW (I) 


Format 

20A4 

615 

5E12.0 

6E12.0 

6E12.0 

6E12.0 

615 

6E12.0 

6E12.0 

6E12.0 

6E12.0 

6E12.0 


* EUher Card 8A £or number densities or Card SB for mole fractions is 
read but no t bo th . 

+ If Card 8A is read, Cards 10 and 11 are not read. If Card 8B is 
read both Cards 10 and 11 are read. 

0 If IEZ = 0 Card 9 is not read. 





TABLE C.2 


VARIABLE DEFINITIONS FOR LRAD 3 


Variable 

TIT 

NETA 

MF 


NS 

LINES 


IDG 


IEZ 


R 

DELTA 

DTIL 

XM0L 


RDZ 


Description 

Title for identification of the problem.' 

The number of points used in the slab calculation. 

Species concentration option variable. 

MF=1 Mole fractions are input on Card 8B and 
number densities are computed. 

MF=0 Specie number densities are input on Card 8A. 

The number of species to be input. 

Line radiation option variable. 

LINES=0 Only a continuum calculation is done. 

=1 A coupled line-continuum calculation is 
done. 

A switch to allow intermediate printout. 

IDG=0 Only finial results are printed. 

=1 Print at each ETA is given. 

-2 Complete print is given. 

The number of points used in the flux integration. 
IEZ=0 The ETA array will be used for the ETZ 

array. 

0<IEZ<NETA Specifies the number of points in the 
ETZ array. Card 9 will be read. 

Body radius (ft.) 

Nondimens ional stand-off distance (6/R) 

Transformed stand-off distance (6) 

A molecular radiation option switch. 

XM0L=O Molecules not included in the radiation 
calculation. 

=1.0 Molecules included in the radiation cal- 
culation. 

The density directly behind the shock. 

(PfijO’ lbn, / ft3 ) 

The temperature profile across the shock layer. (°K) 


T(I) ,1=1, NETA 
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YD (I) 
ETA(I) 
INDEX 


DENS (K 
K= 

FRAC(K 

k=: 

ETZ(I) 

RH0(I) , 
AMW(I), 


I 


TABLE C . 2 (Cont.) 


, I-1,NETA The nondime ns ional shock layer location where tempera- 
ture, concentration etc. are given. (y/6) 

, I=1,NETA The Dorodnitzyn transformed shock layer locations 
corresponding to the y/6 locations. (T)) 

i 

The number given each specie for use in storing arrays. 
This permits species to be read in any order and is 
placed before each set of cards of type 8A or 8B. 


1 = 02 

7 

= H 

2 = N2 

8 

= C2 

3 = 0 

9 

= H2 

4 = N 

10 

= C0 

5 = E- 

11 

= C3 

6 = C 

12 

= C2H 


, INDEX) 
L,NETA 


Species number densities, (partic les/crn ) 


, INDEX) 
. ,NETA 


Species mole fractions. 


1=1, IEZ A subset of ETA points used in the flux integration. 

If IEZ=0 the ETZ points will automatically be set 
equal to the ETA(I) points. 

I=1,NETA Nondimens ional density profile across the shock layer. 

I=1,NETA The average molecular weight profile across the shock 
layer . 
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The species considered in the program will probably be a subset of 
the species considered in a shock layer flow- field. Only those twelve 
species listed with an INDEX number can be input. If the field of 

interest is ionized, the input of electron concentrations is necessary, 
although ionic species are not input. 

The subset of ETA points, the ETZ points, used in the flux 
integration calculation should be carefully chosen. The purpose of 
using a subset of ETA points is reflected in the required computation 
time (i.e. 5 minutes for 60 ETZ = ETA points and 2.5 minutes for 30 ETZ 

points where NETA = 60). The computation times cited are for a IBM 
360-65. To maintain accuracy keep ETA points as ETZ points in regions 
where either concentration or the temperature is varying rapidly. 

OUTPUT DESCRTPTT ON 


This section presents a description of the LRAD 3 program output 
format and definition of output symbols. The reader may find it 

instructive to refer to the listing of the sample problem while 
reading this section. 

The first four pages of output are a print of the input data. 
This is provided so that the user can check the input for possible 
errors, and it also provides identification of the problem. All 
quantities on these pages are defined in the input guide section. The 
concentrations printed on pages 3 and 4 are either number densities or 
fractions depending on the option used for input purposes. The 
difference is easily distinguishable by the magnitude of the numbers. 
These concentrations are followed by the shock thickness DELTA in cm. 
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The standard output shown (i.e. IDG=0) provides a print of the specie 

number densities on pages 5 and 6. The results of the radiation , 

calculation are printed on the seventh page of output. The continuum 

contribution and the line contribution to the spectral flux are printed 

for three ETA points (ETA=0.0=wall, ETA=0.5, ETA= 1.0= shock) as a 

function of frequency interval HNU=hv . The columns of fluxes in 

watts/cm denoted by QPLUS and QMINUS designate fluxes to the surface 

and away from the surface respectively. The total radiative flux at 

2 

the three ETA locations stated above are printed in watts/cm . The 
number on the left can be interpreted as the radiative heating to the 
surface. 

SAMPLE PROBLEM AND PROGRAM LISTING 

The following example is presented to illustrate the basic input 
for the LRAD 3 program and to show a typical output listing. The example 
considered is typical of stagnation line shock layers in that the 
temperature and species compositions change from surface to post shock 
values. All twelve possible species are included in the example. The 
integration is carried out over 59 ETA=ETZ points and thus card-type 
9 is not input. Furthermore, species number densities are input on 
card type 8A and thus card type 10 and 11 are not required in addition 
to RDZ on card type 3. 

The required input data for this example are listed on the 
following pages. After the input a listing of the output for the example 
is given. Finally a listing of the LRAD 3 program is provided. 
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APPENDIX D 


VISRAD 3 'COMPUTER PROGRAM 

i 

DISCUSSION OF THE PROGRAM 

This appendix describes a digital computer program which is 
designed for prediction of stagnation line viscous and radiative 
coupled shock layer structure and the resulting heating rates pro- 
duced by a blunt body during super-orbital entry into planetary 
atmospheres. Formulation of the problem was presented in Chapter 2 
from which the first order stagnation line shock layer equations 
were derived. The computer program described herein was developed 
using the stagnation line equations and the implicit finite difference 
techniques stated in Chapter 4. The regime of atmospheric flight 
is restricted to the laminar continuum regime and to conditions where 
thermodynamic equilibrium can be applied. The governing equations 
include the effects of mass injection of ablation species, radiation 
cooling, and coupling between convection and radiation. The radiative 
transport term is evaluated using either of two models. The first 
one is the emission model described in Chapter 3. The second is a 
line and continuum band model described in Chapter 3 and Appendix B. 

Solutions to this problem are obtained by numerically solving 
the set of stagnation line thin shock layer equations developed in 
Chapter 2. Chapter 4 presents the finite-difference method used for 
both the momentum and energy equations as well as the convergence 
method and criterion. The entire solution process consists of a 
complex iteration scheme which is started by assuming initial values 
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for the solution profiles and then iterating until convergence is 
obtained. 

The VISRAD 3 program utilizes thermodynamic and transport 
properties of air from Ref. D.l and thermodynamic properties of 
air- ablation product mixtures from Ref. D.2 as discussed in Chapter 
3. The calculation of these properties were programmed as separate 
subroutines such that either or both calculations could be made 
during the solution procedure. Air values for 1^, p and Cp T are used 
throughout the shock layer. An elemental step function is assumed 
for the elemental species solution. The elemental switch from abla- 
tion products to air is made at the stagnation point.- Species 
compositions are then computed using the two methods based on the 
elemental composition, local temperature and pressure. 

The program was designed as a tool for use by thermodynamic 
engineers for thermal environment prediction studies. It provides an 
accurate method for analyzing a variety of atmospheric entry heating 
problems. The following is a summary list of the capabilities included 
in the VISRAD 3 computer program: 
o Stagnation line solutions 

o Coupled convective and radiative energy flux calculation 
o Line and continuum radiation calculation 
o Large or small mass injection rates of ablation species 
o Complete equilibrium chemistry 
In addition to these basic capabilities, the program is designed so 
that other options can be added without altering the basic program. 



A solution to the elemental continuity equations for binary diffusion 
has been successfully added to the basic program by Esch (Ref. D.2). 
Additional modifications which could be incorporated into the program 

I 

are different diffusion, transport and chemical models. 

There are three heating rate options available in the- program. 

o Convective heating only 

o Uncoupled convective and radiative heating 

o Coupled convective and radiative heating. 

PROGRAM PROCEDURES 

The VISRAD 3 computer program was developed using a philosophy 
of minimizing a users effort and maximizing program flexibility and 
adaptability. Accordingly, the basic program logic as shown in 
Fig. D. 1 is quite simple. However, these five basic subprograms are 
supported by 20 subroutines and one function subprogram. Each of 
these peripheral programs perform specific computational functions 
thus enabling direct addition or substitution of other subroutines 
when required. 

In order to minimize input requirements, three techniques were 
used. The first consists of internal initialization of values for 
temperature, density, viscosity and stand-off distance which are 
necessary to start the solution procedure. If better guesses are 
available they may be input as discussed in the next section. The 
second technique involves internal specification in BLOCK DATA of 
problem defining parameters such as the elemental composition at the 
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FIGURE D.l BASIC LOGIC DIAGRAM 
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surface and the thermodynamic curve- fit constants which are changed 
quite infrequently. The third technique consists of internally 
selecting program options if an option variable is left blank on an 
input card. In this procedure the most commonly used options are 

performed when a blank is input. 

The use of the elemental step function for the solution of the 
elemental species equation presented one computational difficulty. 

This type of solution produces a jump in average molecular weight at 
the stagnation point which results in a discontinuity in density. 

To eliminate this problem the average molecular weight of the 
ablation products and that of air was linearly weighted with distance 
in 1 \ near the stagnation point. This is an approximation of the effect 
diffusion has on the average molecular weight. In this manner the 
density profile, which is primarily dependent on the temperature, is 
computed as a smooth function. 

The start up of the VISRAD 3 program can be achieved in a number 
of ways. As stated previously internal guesses are available to begin 
the iteration procedure. Two types of temperature profiles are 
available. One for no mass injection and the other for mass injection. 
These profiles are usually quite satisfactory as initial guesses if 
an emission radiation coupled problem is to be run. However, if a 
line and continuum radiation coupled problem which includes mass 
injection is to be run the internal guess may not be accurate enough. 
Consequently, a start up logic has been developed to improve the 
initial guess by beginning the solution procedure using air properties 



IEM (Energy- Momentum Iteration No. 


Selected iteration 
to add ablation 
products to the 



FIGURE D. 2 VISRAD 3 Start-up Option Logic 
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and the emission model. The start up logic is presented schematically 
in Fig. D.2. The computation procedure is begun with internal 
guesses and procedes for n- 1 iterations as if an air-emission coupled 
problem is being solved. At the nth iteration the program automatically 
introduces ablation products and the detailed radiation model into 
the solution procedure. The n-lth iteration solution profiles are 
used as initial guesses for the "new" problem. In this way reasonably 
accurate initial guesses can be achieved with minimal computer time 
for the specific problem under consideration. The user must select 
the iteration number at which the switch is to be made. For most 
problems the simplified solution is nearly established after five or 
six iterations. The iteration number selected can take on any value 
from one up. Thus this start up option can be eliminated by introducing 
the ablation species and detailed radiation at the first iteration. 

If the iteration number is chosen so large that the switch is never 
made an air, emission-only solution is obtained. Fig. D.2 also shows 
the four iteration numbers where a step- size adjustment is made. Since 
the choice of iteration number for switching ablation and radiation 
models is arbitrary, some numbers will allow two step- size adjustments 
to be made at a single IEM iteration (i.e. 2 or 4). 

The procedure for step- size adjustment along with the convergence 
criterion and procedure are presented in detail in Chapter 4 and thus 
are not discussed here. 

INPUT GUIDE 

All inputs to the VISRAD 3 computer program are read from cards 
supplied by the user; no tapes are required. The basic inputs consist 



of free- stream velocity and density, principle body radius, wall 
temperature and mass injection rate. Additional input parameters 
are needed to determine which program options are desired such as 
coupled or uncoupled solution, emission or line and continuum 
radiation model. Four basic formats are used for input, 15, E12.0, 
E15.8 and A4. Multiple case runs, and hence entire trajectories 
can be processed by placing the input data for each new case behind 
the data for the previous one. The entire run is terminated by an 
end of file. 

Tab. D.l presents the format of the card input and Tab. D.2 
provides a corresponding definition of variables. 
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TABLE D.l 

CARD INPUT FOR VISRAD 3 


Card 

Variables 

Format 

1 

TITLE 

18A4 

2 

KEEP, NETA, IRAD, ITYPE, MAXM, MAXE , 
MAXD, LT, IPHI , FPRCT , TPRCT , IDEBUG 

915, 2E12.0, 2X, 11 

3 

UINF, RINF, R, TWK, HT0TAL, RVW 

6EI2.0 

4 

DELTA, DHL , RZB , RE, PDTIL 

5E12.0 

5 

T(I) 

6E12.0 

6A 

RH0(I) 

6E12.0 

6B 

RM(I) 

6E12.0 

7 

DEPS 

E12.0 

8 

ETA(I) 

6E12.0 

9 

NDEBUG , IRS, IAB 

315 

10 

CWALL(J) 

5E15.8 



TABLE D. 2 


VARIABLE DEFINITIONS FOR VISRAD 3 


Variable 

TITLE 

KEEP 


NETA 


IRAD 


ITYPE 


MAXM 


MAXE 


MAXD 


LT 


Description 

Title for identification of the problem. 

Indicator to determine if the temperature 
profile from the previous case is to be kept 
as a guess for the current case. 

KEEP = 0 Temperature not kept 
KEEP = 1 Temperature kept 

The number of points to be used in the shock 
layer profile. If NETA - 0, a set of 51 equally 
spaced points will be used. If NETA > 0 card 8 
must be read. 

A variable used to specify the type of solution. 
IRAD = 1 Convective solution only 

= 2 Uncoupled radiation solution 
= 3 Coupled radiation solution 

A variable used to specify the type of radiation 
model to be used. 

ITYPE = 0 Line and continuum radiation model 
= 1 Emission radiation model 

Maximum number of iterations allowed in the 
internal momentum loop. If MAXM =0, it is 
internally set = 15. 

Maximum number of iterations allowed in the 
energy equation and in the overall momentum- 
energy loop. If MAXE - 0, it is internally 
set = 15. 

Maximum number of iterations allowed in the 
external momentum loop. If MAXD =0, it is 
internally set = 15. 

Indicator to determine if a temperature guess 
and if p & pp. guesses are to be read in. 

LT = 0 Cards 5 and 6 are not read. 

= 1 Card 5 but not card 6 is read. 

= 2 Cards 5 and 6 are read. 



TABLE D.2 (Cont.) 


I PHI 


FPRCT 


TPRCT 


IDEBUG 


UNIF 

RINF 

R 

TWK 

HTOTAL 

RVW 

DELTA 

DTIL 

RZB 

RE 


Indicator to determine if the shock curvature 
is to be input. 

I PHI = 0 de/d§ = 0 is internally set. 

= 1 Card 7 is required for input. 

Convergence tolerance for each point in the t' 
profile. If FPRCT =0.0 it is internally 
set = .005. 

Convergence tolerance for each point in the T 
profile. If TPRCT = 0.0 it is internally 
set = .005. 

A switch to allow intermediate printout to 
be obtained at each iteration 
IDEBUG = 0 No print. 

= 1 Print is given. 

The freestream flight velocity (U^) in feet/sec. 

3 

The freestream density (p^) i- n slugs/ft . 

Principal body radius in feet. 

Wall temperature in degrees Kelvin. 

2 2 

Total freestream enthalpy in ft / sec . If 
HTOTAL = 0.0, it is set to U^/2. (Freestream 
static enthalpy is assumed negligible). 

Mass injection rate (pv) /( pU) . 

W CO 

An initial guess for the shock standoff distance 
6/R. If DELTA = 0.0, a guess is supplied by 
program. 

A guess for the transformed standoff distance 
*6/R. The program will also supply this value 
if DTIL =0.0 

The density ratio across the shock p = p m /p§* If 
RZB is input as 0.0, the code will determine a 
value. 

The Reynolds number for the problem, Re g = 

U Rp./p,.. This quantity is determined by the 
program if RE is input as 0.0. 
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TABLE D. 2 (Cont.) 


PDTXL 


T(I),I=1,NETA 


RH0(I),I=1,NETA 


RM(I) ,1=1 , NET A 


DEPS 


eta(i),i=i,neta 


NDEBUG 


Convergence tolerance placed on 6 for 
solution convergence. If PDTIL - 0.0, 

■ 1 V C* O t* — .001 > 


total 
it is 


An initial guess for the dimensionless shock 
layer temperature profile (T/T^) . If LT>0, 
this profile is supplied by the user. 


An initial guess for the dimensionless shock 

layer density profile ( p /p 5 ) - If Li * 
this profile is supplied by the user. 


An initial guess for the dimensionless 
layer pp profile (pp/p^Pc)- If ^T 
profile is supplied by the user. 


shock 

this 


The stagnation line shock curvature (dc/d§). 
If IPHI = 0 then de/d§ = 0.0 is interna y 
set. If IPHI = 1» card 7 is read and de/d^ 
is supplied by the user. 


The grid shock layer points at whicMhesi 
tion profiles are to be computed. If NETA 0, 
ati is set to 0.02 and ETA(I) is computed by the 
program. (ETA(l)=0.(Hwall, ETA(NETA)=1 .0* shock) 

Debug option to output thermodynamic curve- 
fit equations. 

NDEBUG = 0 No output. 

= 1 Output given. 


IRS 


LAB 


Indicator to determine if the wall species 
composition is to be input on Card 10. 

IRS = 0 Card 10 not read and wall mass 

fractions are set equal to those for 
a nylon- phenolic ablator, 

IRS > 1 Card 10 is read. Wall mass factions 
are provided by the user. 

Indicator to determine at what iteration (IEM) 
the ablation products are included in the 

computation. t , 

IAB = 0 Ablation products are introduced 

at IEM = 5. 

IAB > 0 Ablation products are introduced 
at IEM = IAB - 1. 
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TABLE D. 2 (Cont.) 

CWALL(J) , J=1 ,NSP Wall mass fractions. 


NSP 

= 20 

4- 





J = 

1 ** °2 

6 -* N 

11 -+ 

CO 

16 -4 

C 3 H 


2 -» N £ 

7 -♦ E" 

12 

C 3 

17 -> 

C.H 

4 


3 -♦ 0 

8 -♦ C 

13 -♦ 

CN 

18 -♦ 

HCN 


4 -> N 

9 -♦ H 

14 -♦ 

pa 

CM 

19 -» 

C 2 


5 -» 0 + 

1° -h 2 

15 -> 

C 2 H 2 

20 -» 

C* 
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Some caution should be exercised when preparing an input for this 
program. The present program considers 5 elements, including electrons, 
and twenty species listed under CWALL(J). The set of species was 
selected for an air atmosphere and an phenolic- nylon ablator. If 
another ablator is selected for study and this set of species is 
appropriate, no alteration of the program is required. All that is 
required is a card input of wall mass fractions of the ablator selected 
on Card 10. If extensive study of a different ablator using this 
program is anticipated, the user may find it convenient to change the 
wall composition stated in BL0CK DATA under CWALL rather than reading 
in the data for each run. If required, a change to another set of 
species for thermodynamic calculations can be made with compar itively 
little difficulty. This may be achieved by changing the thermodynamic 
curve-fit constants in BL0CK DATA. The thermodynamic curve- fit 
equations were listed in Table 3.2 and the correspondence between the 
coefficients in the table and those in the program is 
For 



1000 < T < 6000 



6000 < T°K 


AI 

A 1 

ss 

All 


BI 

A 2 

« 

BII 


Cl 

A 3 

= 

CII 

» 

DI 

A 4 

= 

DII 


El 

A 5 

- 

Eli 

• 

FI 

A 6 

= 

FII 


GI 


GII 



where the coefficients are dimensioned to include a value for each 
species. The species ordering is given in the SP array with 
corresponding ordering in SMW (i.e. species molecular weight) array 
in BL0CK DATA. Curve- fit constants in the proper format for over one- 
hundred species are available in Ref. D.3. If a different chemical 
system is chosen for study a check should be made to determine if the 
species added are radiatively important. This may be achieved using 

the method suggested in Chapter 3. 

An accurate initial guess appears to be necessary for some 
flight conditions and body radii to avoid convergence difficulties. 

The most sensitive profile appears to be the shock layer temperature. 

A library of starting values for the temperature and T| profiles should 
be developed. The VISRAD 3 program provides punched output of these 
three profiles along with the eta, 1), step-size profile. The punched 
output is in the format for input stated for card type 5, 6A, 6B and 
8. This not only provides the user with a means of collecting 
guesses but also provides the program with a restart capability. 

An additional input hint which has proven useful is to run using 
the emission model first and use this output as input for a run with 
the detailed radiation model. It has been found that coupled 
solutions using the emission model are not as sensitive to the input 
guess as solutions using the detailed radiation model. This method 
or the start up procedure discussed in the previous section is to be 
used if a library of guesses has not been developed. 
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A similar technique has proven useful in running entire trajectories 
(usually without mass injection). By inputing KEEP = 1, the program 
will use the solution at the Nth trajectory point as an initial guess 
for the N+lst trajectory point solution. In the early portion of a 
typical trajectory the radiative coupling is weak and hence the 
solution is not as sensitive to the initial guess. Since radiative 
coupling varies smoothly over the trajectory, the solution at the Nth 
point provides a reasonable guess for the solution at the N+lst point. 

Two changes in the computer program are recommended if extensive 
use of the program is anticipated. First, substitute a calculation 
of the total heat capacity of ablation products for that of air which 
is presently used. This calculation should be made in subroutine 
PROPRT where the frozen part of the heat capacity for the ablation 
products is presently computed. The calculation for the reacting 
part may be computed using the method of Ref. D.2. Secondly, 
additional logic should be built into subroutine STPSZE which internally 
adjusts the T] step- size. Some numerical difficulties are encountered 
if the step size is increased to AT) = .04 or. the shock side of the 
stagnation point for (pv) w > .10 even though the temperature gradient 
is small. This subroutine should be modified to provide a smaller 
step size in this region. 

OUTPUT DESCRIPTION 

This section presents a description of the program output format 
and definitions of the output symbols. The reader may find it 
instructive to refer to the listing of the sample problem presented 
in the next section while reading this section. 
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The first page of output is a print of the input data. This is 
provided for a check of the input and an identification of the 
problem. All quantities on this page are defined in the INPUT GUIDE 
section. The second page also contains problem specification data 
which is self explanatory. Following the guessed nondimensional 
stand-off distance (DELTA) and transformed stand-off distance (DTIL) , 
a listing of the dimensional stand-off distance computed at each 
iteration is given if the detailed radiation model is used. 

Species number densities for those species used in the radiation 
calculation during the final iteration are printed on the third page. 
The fourth page provides an output of radiative fluxes computed 
during the final iteration. The continuum contribution and line 
contribution to the spectral flux is printed for three ETA points 
(ETA = 0.0 = wall, ETA = stagnation point, ETA = 1.0 = shock) as a 
function of frequency intervals and frequency centers respectively. 

The columns of fluxes in watts/cm^ denoted by QPLUS and QMINUS 
designates fluxes toward the surface and away from the surface 
respective ly. 

The fifth page begins with a cheerful message noting the solution 
converged. Following this message is a printout of the shock stand- 
off distance parameters (DELTA, DTIL) ; the convective (QC) , radiative 
(QR) , and total heating rate; and the radiative heat transfer 
coefficient (CUR). To the left of the heating rate data the density 
ratio across the shock (RB) and the mass injection rate (RVW) are 
stated. Following the heating rate data is a print of the numerical 



441 


of the solution profiles as a function of the shock layer coordinates 
ETA and (Y/D) . The solution profiles printed are: 

F / = velocity function 

RV = pv/o U (nondimens ional mass flux per unit area) 

r r co CO 

T/TD = T/T. (nondimens ional temperature) 

6 

E = (E (radiative flux divergence) 

V = v/U (nondimens ional normal velocity) 

00 

V(fT/SEC) = dimensional normal velocity 

G = nondimensional total enthalpy 

H(STATIC) = nondimensional static enthalpy 

The shock layer thermodynamic and transport properties are out- 
put beginning on the following page. These properties are then 
followed by a print of the species mass fractions output as a 
function of ETA. In addition to the species on the last page of 
output, the mixture heat capacity (CP) and average molecular weight 
(AMW) are also printed. 

The above description is for a standard output. If, however, the 
intermediate print option is used (i.e. IDEBUG > 0) additional infor- 
mation is printed related to the momentum and energy equation 
convergence as well as the information just described given at each 
energy- momentum iteration. The additional information is clearly 
labled and thus no difficulties should be encountered in its 
interpretation. 

SAMPLE PROBLEM AND PROGRAM LISTING 

The following example is presented to illustrate the basic input 
for the VISRAD 3 program and to show a typical output listing. The 
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conditions defining the problem are: 

U = 50000. ft/sec 
00 

p = 2.69 x 10 7 slug/ft^ 

r oo 

R = 9 ft 


(pv) w = .20 

(de/d§)^_ 0 = 0.0 (internally set) 
T = 3450°K 


(internally set) 


w 

C. = .7303 carbon 

1W 

.0729 hydrogen 
.0496 nitrogen 
.1472 oxygen 

Coupled radiative heating rates are to be obtained from the stagnation 
line solution using the line and continuum radiation model. The 
solution is started using a set of T/T^, p/p^, (PM-)/(pM-)g and T| 
profiles and thus LT = 2 on card type 2. Card type 4 is a blank card. 
A listing of the necessary card input data is shown on the following 


page. 
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The output listing for the stated example is listed on the following 

nine pages. A description of the output format is given in the 

previous section. Following the output listing of the sample 

problem a listing of the VISRAD 3 program is given. 

Subsequent to the publication of the following it was noted that 

after card ENER 880 a card stating 
C0E = 2.1 DTIL 

and card ENER 890 should then read 

A1 = -C0EF*(TARM1 + TARM2 - C0E*TARM3) 


9 
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IF(RVW.LE.C.C) GO TC 93 ENEP 1460 
CALL CHEMEC(l.NF) • ENLP 1470 
CALL PRCPRT (NSP* 1 *NF) EMEP 14 M 0 
CONTINUE ENLK 1490 
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APPENDIX E 


RADCOR COMPUTER PROGRAM 


DISCUSSION OF THE PROGRAM 

The computer progrem RADCOR can be used to compute radiative 
heating rates at the stagnation line or around the body. Heating rate 
calculations are based on a radiative cooling parameter correlation 
which is used in conjunction with an isothermal slab radiation 
calculation made in the program. The isothermal radiation calculation 
is made for post shock thermodynamic and species levels. Each point 
on the body is treated independently and thus the heating at the 
point in question is only dependent on the local shock shape, stand-off 

distance and free stream conditions. 

This program provides the capability of rapidly estimating 
radiative heating rates at the stagnation line or around the body for 
no mass injection. Although its primary use is for earth atmospheric 
entry, heating due to entry into Mars or Venus atmospheres may also be 
computed. The effects of ablation products may be accounted for by 
hand calculation methods described in Chapter 5 and 6. Thus the 
philosophy in using this program is one of obtaining preliminary design 
estimates of the radiative heating environment about a vehicle. 

The basic assumptions of the calculation are. 

1. The shock layer can be approximated locally as an infinite 
plane slab. 

2. The body is assumed spherical. 

3. The shock wave is assumed concentric. 
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4* The shock stand-off-distance is computed using 

6 = p/ (V 8 P/ 3 + ^ 

5. Line and continuum radiation of species H, C, 0 and N are 
inc luded. 

6. Continuum radiation of species CO, 0 2> C 2> N 2> C 3 and H 2 
are included. 

7. Radiation blockage by ablation species is not included. 

8. The surface radiative flux can be computed from the isothermal 
flux using a radiative cooling parameter correlation. 

The principle option of the program computes the surface radiative 

heating for an air atmosphere using the following computation sequents. 

First, the Rankine-Hogoniot equations (2.79 to 2.82) are solved using 

the air thermodynamic properties of Hansen (Ref. E.l). Second, the 

shock stand-off-distance is computed. Third, using the stand-off- 

distance and post shock species composition from the solution of the 

Rankine-Hugoniot equations, the isothermal radiative flux is computed. 

Fourth, the radiative cooling parameter, T , and surface radiative 

heating, q_ , are computed using the following relations from Ref. E.2. 
R 

(q R ) 

isothermal .04 < T < 1.0 (E-D 

X - 3 

% p u 

2 K c o CO 


q R = (0.2 - 0.295 log 10 D(q R ) 


(E.2) 


isothermal 

If the around the body option is used this sequence is repeated at 

two degree increments around the body. 

If the non-air atmosphere option is used the density ratio across 
the shock and the post shock species, temperature pressure and average 
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molecular weight must be input and only stagnation line calculations 
can be made. This option is appropriate for bodies entering either a 
Mars or Venus atmosphere. The only apparent deficiency for this type 
of calculation is the possible radiative effect of CN which may 
exist in significant quantities for some flight conditions. 

There are three features of the computer program that make it 
attractive for preliminary estimates of surface heating rates: 

1. The only inputs that are required for air atmosphere problems 
are the free- stream flight conditions and body radius. 

2. A large number of solutions can be obtained in a small amount 
of computer time (less than .10 min of IBM 360-65 time per 

case) 

3. The program has a modest computer storage required ( 19 K 
words). 

The next section presents the details of the input procedure. 

The usefulness of this program would be enhanced if a correlation 
for the effects of ablation product injection were included in the 
program. An attempt was made without success to find a suitable 
correlation as discussed in Chapter 5. If however a suitable corre 
lation is obtained in future work, its addition to this program would 
improve present accuracy of making rapid preliminary design calculations 

INPUT GUIDE 

All inputs to the RADCOR computer program are read from cards. 

The basic input consist of the free- stream velocity , free stream 
density _or post shock pressure, and the body radius. Four inputs 
formats are used for input, 15, E12.0, E10.0, and A4. The floating 
point formats must have the decimal point punched on the card. 
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Multiple cases for the same kind of problem may be run. The first 
card input signals the kind of cases to be run. For subsequent 
cases card 1 is not read and the cases are to be subsequently stacked. 
Card 1 of the input also serves to indicate multiple runs of a 
particular kind. Parametric studies may be made by changing the body 
radius or free- stream velocity holding the other specifying variable 
constant. This is accomplished by setting NR = no. of body radii or 
NV = no. of free- stream velocities to a value greater than 1. If this 
is done, the base case specified by either card 3A or 3B is computed 
first then the velocity and or body radius is sequentially increased by 
constant intervals according to the following logic. 


R = R. .„. , + 2.0 K (ft) 
initial 


where K = 2 , . ...» NR 


U = U + 2000. J (ft/sec) 

°° “initial 

where J ~ 2, ...., NV 

The U updating takes place in a loop internal to a loop updating R. 

CO 

Consequently for a given R all NV velocity cases will be run. If the 
R is updated the velocity is reinitialized and the sequence is repeated. 

Table E.l provides the details of the card input and Tab. E.2 
provides a corresponding definition of variables. 
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TABLE E.l 

CARD INPUT FOR RADCOR 


Variables 

ITYPE, IB0DY, NR, NV, IATM 
TITLE 

UINF, PDK, R 
UINF, RINF, R 

NETA, LINES, IDG, IEZ, XM0L 
TD, AMW(NETA, RZB , RE, PD 
FRAC(NETA,I) 


Format 


18A4 

3E12.0 

3E12.0 


6E12.0 

6E12.0 


* The value of ITYPE determines which card rs to be read. 

ITYPE = 0 Card 3B is read 
. 1 Card 3A is read 

$ If this card is blank the variable are internally set. 

NETA = 5 
LINES = 1 
IDG = 0 
IEZ = 5 
XM0L = 1.0 

# The value of IATM determines whether these two cards 

IATM = 0 Neither card is read 
= 1 Both cards are read 
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Var i ab le 
ITYPE 


IB0DY 


NR 

NV 


IATM 


TITLE 

• UINF 


TABLE E. 2 

VARIABLE DEFINITIONS FOR RADCOR 


Description 

Indicator used for the specification of which option 
is used for normal shock calculations. 

ITYPE = 0 Shock conditions computed from free 

stream velocity and free stream density . 

= 1 Shock conditions computed from free 

stream velocity and post shock pressure . 

Counter used for the- specification of whether solution 
is found at the stagnation line only or for an around 

the body concentric shock. 

IB0DY = 0 or 1 Stagnation line solution only. 

> ^ Around the body solutions for a 

concentric shock. Value of IB0DY 
determines how far around the body 
the solution is carried. 

Counter determining whether R is to be internally 
updated. 

NR = 0 R is read for each case. 

si R is to updated internally. Value of NR 

determines the number of times R is updated. 

Counter determining whether UINF is to updated 
internally 

NV = 0 UINF is read for each case. ^ 

^ 1 UINF is to be updated internally. Value of 
NV determines the number of times UINF is 
updated . 


Counter determining whether solution is for air 
atomosphere or any atmosphere. 

IATM = 0 Solution for air atmosphere only. 

= 1 Solution for any atmosphere. This 

solution requires the inputing of the 
mole fractions for species present in 
the system. 

Title for identification of problem. 

Free stream velocity (U^ - ft/sec) 

Post shock pressure (P & - atm.) 


POK 
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RINF 

3 

Free stream density (p^ - slugs/ft ) 

NETA 

The number of points used in the slab radiation 
calculation 

LINES 

Line radiation option variable 

LINES 0 Continuum calculation only. 

=1 Coupled line and continuum calculation. 

IDG 

Switch to allow intermediate printout 
IDG - 0 Only final results printed. 
. 1 Print at each ETA is given. 
-- 2 Complete print is given 

IEZ 

The number of points used in the flux integration. 

IEZ = 0 The ETA array will be used for the ETZ array 
0 <IEZ<NETA Specifies the number of points in the 
ETZ array. 

XM0L 

A molecular radiation option switch 

XM0L = 10" ^ Molecules not included in the 
radiation calculated 

= 1.0 Molecules included in the radiation 

calculation. 

TD 

Post shock temperature (T^ - K) 

AMW(NETA) 

Average molecular weight of all species present in 
the system being analyzed. 

p 

RZB 

t ■\ 

Density ratio across the shock Cl. 

RE 

Reynold's number = Re & = 

PD 

Post shock pressure (P § ) 

FRAC(NETA,I) 

Post shock mole fractions of species I used in 
input of non-air atmospheres. 

1 = l = 0 2 4=0 7 = H 10 = CO 

2 = N 2 5 = E" 8 = C 2 11 = C 3 

3 = 0 6 = C 9 = H 2 
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OUTPUT DESCRIPTION 

This section presents a description of the output forest and 
definition of output symbols. The reader may find it instructive to 
refer to the listing of the sample problem given in the next section. 

The first page of output consists of a print of the title card 

•f-iori nn card 4 This provies an identification 

and the options specified on card 4. in f 

of the problee and a check on soee of the pertinent options. 

The second page of standard output (i.e. IDG - 0) is a print of 

the results for one case. Under the heating of "SHOCK LAYER CAS 

PROPERTIES" the free-streae and post shock conditions are given. Names 
and eeanings not listed in the foregoing section are: 

HT0TAL = Free- stream total enthalpy 
VD = Post shock normal velocity 

(R*U)INF = Free- stream mass flux per unit area 
Under the heading of "SPECIES HASS FRACTIONS" the post shock species 
mass fractions of air are listed if the air atmosphere option is used. 
No listing is given for other atmospheres. Under the heading of 
"RADIATIVE FLUX PROPERTIES" information pertaining to the radiative 
heating calculation is given. The names and descriptions of the 
variables printed are: 


PATHLENGTH 

GAMMA 

ISOTHERMAL FLUX 
ACTUAL FLUX 

CHR 


= Path length used in the isothermal radiation 
calculation which equals the stand-off 
distance - 

* P of Eq. (E. 1) 

_ ) of Eq. (E. 1) . 

ISOTHERMAL 

= (surface heating rate) 

= Radiative heat transfer coefficient 
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DELTA/R = 6 (nondimens ional stand-off distance) 

3 

„5RINF*UINF3 = .5 p U (free stream kinetic energy flux 

00 CO 

per unit area) 

qrr «* (qT,)/(q D )n the heating rate referenced to 

the stagnation value 

If the parameter IDG is greater than 0 additional output from 
the radiation calculation is given. This output is similar to that 
discussed in Appendix C and thus will not be detailed here. 

SAMPLE PROBLEM AND PROGRAM LISTING 

The following example is presented to illustrate the basic input 
for the RADC0R program and to show a typical output listing. The 
conditions defining the problem are 
U = 50000 ft/sec 

CO 

p = 9.0 x 10 7 slug/ft"^ 

CO 

R = 9 ft 

The heating rate results shown are for the stagnation point in an air 
atmosphere. The following is a typed listing of the necessary input 
cards for this problem (the zeros need not be punched). 


0 0 

0 

0 

0 

SAMPLE CASE 




50000. 

9.0E7 

9.0 


0 0 

0 

0 

0.0 


The four cards required for this sample case are cards 1, 2, 3B and 4. 
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A computer output listing for this problem is given on the 
following pages. This output listing is followed b y a Fortran listing 

of the computer program. 
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SHOCK LAYER GAS PROPERTIES- 
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